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ABSTRACT 

The piecewise parabolic method and related schemes are widely used to model 
stellar flows. Several different methods for extending the validity of these meth¬ 
ods to a general equation of state have been proposed over time, but direct 
comparisons amongst one-another and exact solutions with stellar equations of 
state are not widely available. We introduce some simple test problems with 
exact solutions run with a popular stellar equation of state and test how two 
existing codes with different approaches to incorporating general gases perform. 

The source code for generating the exact solutions is made available. 


1. Introduction 


Many stellar flows are modele d wi th com pressib le h ydrod ynamics methods, with the 
piecewise parabolic method (PPM) flColella fe Woodwardlll9841) being one of the most pop¬ 
ular algorithms. All extant astrophysical hydrodynamics codes that include a general equa¬ 
tion of state augment the PPM algorithm with additional information about the energy of 
the fluid and use Riemann solvers designed to approximate the wave structure of the solution 
without requiring expensive equation of state calls. However, little attention in the literature 
has been given to verifying these simulation codes directly for the degenerate stellar equa¬ 
tion of s tate. While plenty of shock-tube Riemann problems exist for a simple gamma-law 
gas (e.g. ISodl 119781 ). these do not test the algorithms when the gamma of the gas changes 


dramatically in the problem. We introduce some shock tube problems for the general stellar 
equation of state and compare the hydrodynamic solution to the exact solution of the Rie¬ 
mann problem for these tests. Different codes make different choices when dealing with a 
general equation of state, but no comparisons and little verification has been done on these 
implementations. The basic verification of these methods for stellar flows is the goal of this 
paper. 
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An important diagnostic quantity is the temperature of the gas in the star. Degen¬ 
erate gases are weakly sensitive to temperature—this creates interesting challenges for hy¬ 
drodynamics methods. Traditionally, pressure, not temperature, plays a central role in the 
construction of the fluxes through the interfaces in Godunov schemes for com pressi ble hydro¬ 
dynamics. In the piecewise parabolic method (PPM) implementation (IColella fe Woodward 
1984T ). the primitive variables, p,u,p are reconstructed as parabolas in each zone and we 
then trace under these profiles to find the information that can reach the interface over 
the timestep. Any small errors in the thermodynamic consistency of the interface state, 
introduced from the extrapolation and limiting procedures used in these methods, can lead 
to a large error in the temperature calculated from the equation of state for a degenerate 
gas. While not needed in pure hydrodynamic flows, temperature is an essential quantity for 
reacting flows and radiation hydrodynamics. For this reason, we look at the performance of 
the different methods in preventing artihcal undershoots and overshoots in temperature. 


In the tests below, we use the publicly-available CASTRO code (lAlmgren et ahl 20101 ) 
as our reference and implement the variations as options in CAS TRO. We also compare to 
the original dimensionally-split PPM s olver in FLASH (Fr vxel l et ah 20001) . that in turn was 
based on the PPM implementation in Frv xell et ah (119891) . While we note that FLASH has 
other solvers available, this is the default solver in FLASH 4.0. These co des use s lightly differ¬ 
ent im plementations of PPM, based on original m ethod as describe d bvIColella fe Woodward 
( 1984 ) (henceforth CW) and the description in Miller fe Colella ( 2002 ) (henceforth MC). 
They also take different approaches to the extension of the algorithm to a general equation 
of state. Finally, we note a small correction to the scheme used in FLASH, in Appendix I A. 3 1 


2. PPM Overview 

The conservative Euler equations appear, in one dimension, as: 



dp d(pu) 
dt dx 

= 0 

(la) 

d(pu) 

dt 

d(pu 2 ) dp 
dx ^ dx 

= 0 

(lb) 

9(pE) 

dt 

d(puE + up) 
dx 

= 0 

(lc) 


where p is the density, u is the velocity, p is the pressure, and E is the total specific energy. 
The equations are closed via an equation of state, 


p = p(p, e ) 


( 2 ) 
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where e is the specific internal energy, e = E — \u 2 . For many astrophysical equations of 
state, temperature (T), not energy, is an independent variable, and we have 

P = P(P,T), e = e(p, T) (3) 


and iterations (usually via the Newton-Raphson method) are done to find the T that gives 
the desired e, allowing us to then find p. 

The hydrodynamics equations in terms of the standard set of primitive variables, p, w,p, 

are: 


dp dp du 

m +u ai + p di 

= 0 

(4a) 

du du 1 dp 
dt U dx p dx 

= 0 

(4b) 

dp dp 2 du 

di +u di +pc di 

= 0 

(4c) 


where c is the speed of sound defined in terms of the adiabatic index Ti = <9(logp)/<9(logp)| s 
as 

c 2 = T lP /p (5) 

As most papers deal with a constant-y ideal gas equation of state, it is important to note that 
the pressure equation above is completely general, and can be derived by writing p = p(p, s), 
and differentiating along streamlines with no entropy sources. No assumptions about the 
constancy of 7 are needed. 


For a general EOS, we need to augment our system with additional thermodynamic 
information for the Riemann solver. CASTRO adds the evolution of (pe) to the system: 


d(pe) , d (pe) 
dt +U dx 


, du 
ph Tx 


0 


( 6 ) 


where h = e + p/p is the specific enthalpy. Alternately, Colella fc Glazl ( 19851) (henceforth 
CG) define y e = p/(pe) + 1 and derive an evolution equation for it. Differentiating, we find: 


£le = D (P A _ P D(pe) 1 Dp 
Dt Dt \pe ) {pe) 2 Dt pe Dt 

= (7=-l)(7.-ri)^ (7) 

where we used Eqs. fl4cl [5j [ 6 ]) to simplify. If we eliminate du/dx in favor of Dp/Dt using 
Eq. |4cl then we arrive at the expression from CG. We see that for an ideal gas, since 7 = r 1; 
we have D'y/Dt = 0. 
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An analysis of this system (see any standard hydrodynamics text, e.g. Toro 119971 ). 


including one of the energy equations, shows that there are three characteristic waves that 
carry jumps in the variables at the speeds: A-~) = u — c, = u, and A*- + ^ = u + c (with the 
addition of the (pe) or y e equation, the u eigenvalue becomes degenerate). We will use these 
superscripts, (—), (o), and (+), to distinguish between the waves in the following discussion. 


The PPM algorithm advances the solution through a timestep using a finite-volume 
framework. In each zone, the average value of the state is changed by constructing the fluxes 
through the interfaces of the zones. The update appears as: 


1. Construct limited parabolic profiles of the primitive variables, q l: in zone i using the 
procedure from CW. Thi s result s in a limited parabola, q(x). We note that alternate 
limiters exist (jColella fe Sekoral 1200a ) and are implemented in CASTRO, but we do 
not consider those here. 


2. Integrate under each parabola to find the average state that can be carried to each 
interface by each of the characteristic waves. We define an integral 2+ +;) which is the 
average of the parabola profile of q { to the left (‘—’ subscript) or right (‘+’ subscript) 
of the cell for the region swept out the wave ( v ) (see Figured])- These integrals are: 

(8a) 
(8b) 

with <jC) = \X^\At/Ax (see MC for a discussion and motivation). 


i bta) = 
zPfe) = 


r i +1/2 


crO Arc 

1 


Vi+ i/2-<r ( " ) Ai 


q(x)dx 


a 




q(x)dx 


Vi-l/2 


3. Perform characteristic tracing to build the interface states. Following CW and MC, 
the data in zone i can be used to build the states on the right interface, q^lj^Li anc ^ 

on the left interface, i, q^xj^Ri °f zone i (see Figure [2] for an illustration of the location 
of each of these states). These appear as: 


71+ 1/2 

di+l/2,L 


n+ 1/2 


i+- V C- 

y;A<»>0 

+ 

I 

\ r? ] 

(9a) 

i- A C- 

(?- 

r i U) 

(9b) 


This shows that the same set of eigenvectors and eigenvalues are used for these two 
states. 
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Fig. 1.— Tracing under the parabolic profile to construct l+\q) at the i + 1/2 interface. 


Here, q + and g_ are the reference states, designed to minimize the amount of work 
done by the characteristic tracing (see the discussion in MC). The key part of this step 
is that the sum (Eqs. [9a] and l9bl) only includes the contributions for waves that are 
moving toward the interface. If all waves were included, then this would be a no-op 
(although, see the discussion in [Stone et ah 2008 for HLL-type solvers). 


There is a choice in the primitive variable system one can use, and therefore, the 
resulting left and right eigenvectors: if''* and rf'*. This comes into play with a general 
gas since FLASH and CASTRO use a different “energy” variable in the eigensystem. 
Appendix [A] discusses the options and gives the form of the eigenvectors, and writes 
the states qf+f/ 2 { L r} f° r each variables in a notation similar to that in CW. 


4. Solve the Riemann problem to find the fluxes through each interface and convert the 
unique interface states back to the conserved variables. There are several different 
Riemann solvers used for non-ideal gases in the literature. We discuss some of the 
options later. 


5. Update the conservative variables using these fluxes. 


As a side note, we remark that much of the PPM method (in one-dimension) is done 
to third order, but the conversion from conservative to primitive variables is only second- 
order accurate, and t his c onver sion itself can introduce some thermodynamic inconsistency. 


McCorauoda 

e & Colclla 

(2011 

) discuss how to derive a 4th-order accurate method. Finally, 

we note that 

Colclla 

(19901 and MC described the extension to multiple dimensions. 
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Fig. 2.— The two interfaces states constructed from the data in zone i. 


2.1. Variations Among PPM Implementations 


CASTRO was originally written to closely follow the PPM method as described in MC, 
which differs from the original PPM implementation (CW), and the original dimensionally- 


split PPM solver implemented in FLASH in several ways. Since the original paper flAlmgren et al. 


2010 l). a number of small changes were introduced, which are summarized here. 


• CASTRO originally set the reference state, q ±, to be the cell-center value. It has been 
updated to use the prescription from CW, which uses the integral under the parabola 
for the fastest wave (if the wave moves toward the interface) or the limit of the parabolic 
interpolant if the wave is not moving toward the interface. These options are controlled 

in CASTRO through the castro. ppm_ref erence and castro. ppra_ref erence_edge_limit 
runtime parameters. 

• In the characteristic projection, CW use r = 1/p instead of p in the eigensystem (see the 
form of /3° in CW Eq. 3.7). MC use p. This is discussed in Appendix I A. 21 and explored 
in more detail later. Note however that CW do the parabolic reconstruction of p, not r. 
This choice is controlled in CASTRO by the parameter castro. ppra_tau_in_tracing. 

• In the characteristic tracing, numerous terms appear of the form /? = (/• A q), where 
A q is some jump in the primitive variable q. CW evaluate these /3’s and any other 
quantities in the interface state construction using the reference state (this is equivalent 
to evaluating the eigenvectors with the reference state). MC do not explicitly write out 
the expansion of the dot product and leave things in terms of eigenvectors which are 
constructed using the time-level n data. We use the CW method here in CASTRO, 
which can be selected through the castro ,ppm_reference_eigenvectors. 
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All variants described above use flattening to prevent shocks from becoming too thin, 
but there is variation in when this flattening procedure should be applied. In MC, 
the edge states qi and q r are constructed, then the PPM monotonization procedure is 
applied, and finally the flattening is done to each parabola. In the FLASH split PPM 
solver, the flattening is done on the qi and q r before the PPM monotonization. CW 
do not seem to indicate which of these orderings is preferred. Both of these do the 
flattening before the parabolas are integrated. We note that the flattening function in 
CW is written in a different form than in MC, but they are analytically equivalent. 
CASTRO by default applies the flattening after the integrals X are constructed from 
the parabola, when constructing the final edge states. These differences do not seem 
to influence the solution much. The different methods can be explored in CASTRO 
through the parameter castro. ppra_flatten_before_integrals. We n ote t hat not 
all PPM implementations use flattening—by default ENZO OBryan et al. 2(LL4) has it 


disabled. In the tests that follow, we adopt the FLASH ordering of flattening and then 
limiting in CASTRO. 


Finally, CW use a contact steepening algorithm to keep contact discontinuities thin. 
MC and CASTRO do not implement this, but the FLASH split PPM solver does and uses 
it by default. For the FLASH runs presented here, we disable contact steepening. Some 
authors (e.g. [Stone et al.1 120081 ) suggest limiting on the characteristic variables themselves. 


This is the default in the latest version of the FLASH dimensionally-split PPM solver, and 
we will look at its influence when we compare with FLASH. 


2.2. Extension to a General EOS 

CW do not address how to extend PPM to a general equation of state. Most works 
cite CG as the inspiration for dealing with a general EOS, in particular, for the Riemann 
problem, but other variations exist. CG predict an interface value of y e = p/(pe) + 1 using 
an equation that captures the thermodynamic evolution along streamlines (similar to Eq. [TJ). 
We note that CG did not discuss how to use parabolic reconstruction of the fluid quantities 
with a general EOS. Their y e does not participate in the characteristic tracing described 
above; rather, it is constructed from the predicted interface value of p, using the average Ti 
and 7 e on either side of the interface. Since Vi does not explicitly appear in the fluxes, CG 
argue that taking the cell-average value for the interface is second-order accurate. 

The FLASH dimensionally-split PPM solver reconstructs y e and I i as parabolas and 
integrates under the A*' 0 ) wave to get their interface values (note: this does not appear to be 
documented in the FLASH paper). This differs substantially from CG, and suggests that they 








both should obey a hyperbolic PDE. We show in Appendix IA.3I that it is possible to include 
7 e in the characteristic projection and it can actually jump across all three characteristic 
waves, not just the A^°' wave as assumed in FLASH. We know of no method to allow for a 
high-order reconstruction of to the interface. In the construction of the interface states, 
FLASH evaluates the Lagrangian sound speed, C, using p, p, but with the cell-centered Ti. 
This is correct if we choose not to predict Ti to the interfaces using parabolic reconstruction, 
but potentially inconsistent with a general EOS if you do. In CW, Tx = y e = 7 was constant, 
so the issue of what Tx to use in C was not discussed. 


CASTRO includes an evolution equation for (pe), which appears as an additional hyper¬ 
bolic PDE in the primitive variable system. This is reconstructed and enters into the char¬ 
acteristic tracing in the same fashion as all the other primitive variables (see Appendix lA.ip . 

Solving the Riemann problem exactly for a general equation of state is expensive (see 
Appendix [B]). Both the FLASH dimensionally-split PPM solver and CASTRO use an ap¬ 
proximate two-shock Riemann solver. Here, the left and right waves are assumed to be 
shocks and jump conditions are used to link to the state in-between these waves. This mid¬ 
dle state is traditionally called the “star”-region. To describe the thermodynamics of our 
general gas, auxiliary information is needed to supplement the primitive variables. CG use 
the same evolution equation that predicted the interface value of y e to predict the value of 
7 e in the star-region of the Riemann problem. This is an iterative method that converges 
to find a value of the nonlinear wavespeeds across the (presumed) shocks in the Riemann 
problem. This is then used to construct the energy in the star-region. FLASH uses the 
parabolically-traced values of q e as input to the CG Riemann solver. 


CASTRO instead uses the interface values of (pe) along with p, u, and p as inp ut to the 


Riem ann problem. The default solver in CASTRO follows the procedure in flColella et ah 


19971 ; IBell et ah 19891 ) where jump conditions are used to estimate a value of (pe) in the star 
region. However the option exists to use the CG Riemann solver with interface values of q e 
set as: 

leg = /V + 1 ( 10 ) 


(pe) 

where p s and (pe) s are the predicted interface values of p and (pe), and s G {7,r} indicates 
whether we are to the left or right of the interface. 


Th ere are other potential choices for Riemann solvers. The HLL solvers flHarten et a l. 


19831: Einfel dtl 119881 ) use the jump conditions and simple estimates of the wave speeds to 
give the fluxes directly, and should work for a general gas without issue. The HLLC solver 
advocated by Toro ( Torol 119971 ) would need modification for a general gas because of the 
construction of the speed of the contact wave. 
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Each of these methods incorporates information about e in some fashion. We note that 
there is a potential for thermodynamic inconsistency on the interfaces—the quantities p, 
p, and 7 e or (pe) were brought to the interfaces independent of one another, subjected to 
limiting, flattening, and characteristic tracing. In fact, they are not independent, and there 
is an error in how well these interface states obey the equation of state. 


3. Numerical Tests 


Little, if any, comparison of astrophysical hydrodynamics codes to exact solutions of 
the shock tube problem has been done. Generating exact solutions for an arbitrary equa¬ 
tion _of state is strai ghtforward—these are the exact solutions to the Riemann problem. 
Colella fe Glazl (119851 ) outline the procedure to exactly solve the Riemann problem. In Ap¬ 
pendix [B] we summarize some of the implementation details and give the initial conditions 
for four problems: a Sod-like problem, a double rarefaction, a strong shock, and conditions 
mimicking the edge of an under-r esolved star . The first three are our stellar EOS analogs to 
the standard test problems from Toro (119971 ). The tests presented here are not exhaustive, 
but sample some of the flow conditions we might encounter in large-scale simulations. We 
note that y e is not constant in these tests—and in test 4, it varies through its entire valid 
range. 

We use CASTRO as the main code to test our ideas, but we also run the same tests 


“as-is” with the FLASH dimensionally-split PPM solver (iFrvxell et al. 2000 ) (version 4.0) 
just for comparison. This solver in FLASH can be thought of as the CW method, performing 
a reconstruction of both q e and Ti. In the discussion below, we use “FLASH” to mean this 
specific dimensionally-split PPM solver. 


CASTRO is primarily a 2- and 3-D code, so we do all the runs in 2-D, with the transverse 
direction simply replicating the shock tube problem. This means that the transverse flux 
difference will have no contribution to the inter face states . For all CASTRO tests, we use 
the iterative Riemann solver described in IColell a fe Glaz (119851) (this is enabled with the 
CASTRO option castro ,use_colglaz=l) unless otherwise specified. 


The main driver in FLASH is written with dimensional splitting in mind, so the PPM 
solver takes a cycle consisting of two timesteps with equal At (swapping the ordering of 
the directional sweeps in multiple dimensions) and then re-evaluates At based on the flow 
conditions in the domain. We run FLASH in 1-D, and to have the timestepping match 
CASTRO, we allow FLASH to only take a single timestep per cycle. Both codes will modify 
the last timestep to end at the time defined by the problem. This change to the timestepping 
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is the only change we made to FLASH. 


We note that some variation of these results can be expected with varying CFL number 
or using an initially small timestep, but we will not explore these details. All runs use 128 
zones, a CFL number of 0.8, and take an initial timestep of 0.1 of the CFL step. We also 
restrict the maximum increase in At from one step to the next to be 10%. 


_ Bo th codes use t he sam e full stellar equation of state described in iTimmes fe Swestv 

( 2000 ); Frvxell et ah ( 2000 ). A low temperature floor, T small = 10 4 K, is imposed by the 
equation of state. It is important to note how this small temperature floor is applied. 
During the hydrodynamics solve, we typically enter the EOS with p, e want ,Afc and ask for 
p, T. The equation of state most naturally deals with p, T, X k as inputs, so we must do 
a Newton-Raphson iteration to get the energy, e want we want. During the iteration, if 
T < T sma n, then we stop our iteration and return the thermodynamic state corresponding to 
the temperature floor. The problem is that e sma n = e(p, T sma n, X k ) % e want , so we either break 
energy conservation by resetting the hydrodynamic state to e sma n or break thermodynamic 
consistency by keeping e want . For the simulations presented here, we found that keeping the 
thermodynamic consistency is more important to getting a good temperature field, and this 
is the procedure used in both codes (this is the default behavior in FLASH and controlled 
by the eos_input_is_constant parameter in the feextern namelist in CASTRO). 


3.1. Comparison of Riemann Solvers 


Our first comparison is to loo k at the difference in the solution between the CG Riemann 
solver and the solver presented in IColella et al.l (119971) (henceforth called CGF), which was 
the default Riemann solver used in Almgrcn et ah ( 2010 ). We focus on tests 1 and 4, 
since those are the only tests that show substantial differences. Figure [3] shows the density, 
velocity, pressure, and temperature fields for the two cases together with the exact solution. 
Figure [4] shows the error against the exact solution for these two runs. The temperature 
field shows the most differences (amongst the two solvers and against the exact solution). 
We see an oscillation in the CGF temperature solution behind the contact discontinuity 
and a slightly more pronounced undershoot in the temperature in the CG solver. Both are 
unphysical. Looking at the density and pressure errors, we see slightly larger error between 
the contact and rarefaction in the CGF case than in the CG case. Figure 0 shows the 
comparison for test 4. Here we see differences in the velocity on the right of the domain and 
differences in achieving the cool temperature behind the rarefaction. The CGF gets closer 
to the temperature minimum. Going forward, we use the CG solver exclusively to avoid the 
oscillations see in test 1, but these results leave the door open to exploring other Riemann 
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solvers for general equations of state. 


3.2. Using p vs. r in the Characteristic Tracing 

Here we examine the effect of using r = 1/p instead of p in defining the eigensystem 
used for the characteristic tracing. In essence we are comparing the original CW method to 
CASTRO’s implementation of the MC method. Again, only test 1 shows any difference, and 
again it is slight. Figure [6] shows the temperature results for test 1. The difference here is 
so slight that it seems that this variation does not really matter. 


3.3. Comparisons with FLASH 

Finally we summarize by comparing CASTRO’s method to FLASH. In some sense, 
this allows us to investigate the differences between using y e vs. pe in the interface states 
(although, see Appendix IA.3|h We use both the default options in FLASH (characteristic 
limiting) and a run with the characteristic limiting disabled. Figures 171 through [TUI show the 
results. In test 1, the CASTRO solution appears to suffer the least amount of temperature 
undershoot behind the contact, with the FLASH solution without limiting on the charac¬ 
teristic variables also looking good, and the default FLASH having the deepest undershoot. 
This suggests that the characteristic limiting may not always be a good idea, and its use 
should depend on what quantities you are looking to optimize. For tests 2 and 3, there 
do not seem to be any major differences—all solvers show the same features. For test 4, 
all of the methods have difficulty resolving the narrow region between the rarefaction and 
the contact, and largely look the same, except for the velocity, where all the methods have 
trouble finding the correct velocity at the right edge of the domain. 


4. Discussion 

While the basic PPM algorithm is now 30 years old, there are a large number of vari¬ 
ations in the literature. We looked at the implementations in CASTRO and FLASH, and 
performed verification tests on whether they can accurately model stellar flows. The main, 
simple message of this paper is that it is possible to do verification of stellar hydrodynamics 
codes on problems with a real stellar equation of state, and of course, that one should be in 
the habit of doing this. Part of the inspiration for this exploration was the common appear¬ 
ance of temperature “funniness” such as flooring or oscillations (see, e.g. the comparisons 




12 


in 


Almgren et ah 20060 in hydrodynamic flows. These tests can be used to help design new 


methods that optimize for a particular behavior. 


Our tests showed that the different approaches that FLASH and CASTRO take toward 
incorporating the real equation of state into the solution appear reasonable and consistent 
with one another. They also showed that the updates to CASTRO highlighed in § 12.11 
perform well. In developing the structure of the eigensystem when y e was added, we made a 
recommendation for the FLASH dimensionally-split PPM solver to improve the prediction 
of 7 e to the interfaces. We stress that the tests here are not complete and much more 
exploration can be done. 


The multidimensional treatment of these methods, a nd de aling with a general equation 
of state, is more complex. In an unsplit method (e.g. ICol ellal Il990l ) the interface states 
are first predicted in primitive variable form as if the system were one-dimensional, then 
converted back to conservative form where a transverse flux difference is added. This makes 


the interface states 


see 


what is happening in the transverse directions. Finally, the states 


are converted back to primitive form for the Riemann solve. This conversion can affect 
thermodynamic consistency, especially if it is done without involving the EOS. We note that 
with either (pe) or y e as an auxiliary variable, a transverse term will need to be incorporated. 
This is done in CASTRO currently, and an equation of state call can be avoided by using the 
multidimensional pressure equation throughout the procedure, but a detailed exploration of 
alternatives may still be illuminating. 


Finally we suggest that future stellar hydrodynamics codes test against the exact solu¬ 
tion to the shock tube problem for the stellar equation of state. We make our exact Riemann 
solver and all of the CASTRO improvements tested here freely available in the CASTRO 
code distribution (the Castro Sod_stellar/ problem has all the necessary inputs). Similar 
test problems can be designed for the nuclear equation of state used in core-co ll apse super¬ 
novae simulations, although we note that the solution procedure from lColella fe Glazl (119851) 
we follow does not work for a non-convex equation of state, so modifications will be necessary 
for a general nuclear equation of state. 


This research was supported by NSF award AST-1211563. We thank Ann Almgren, 
John Bell, Alan Calder, Sean Couch, and Dongwook Lee for very helpful discussions and 
feedback. 
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A. Eigenvectors and the Characteristic Projection 

Here we summarize the eigensystems corresponding to different sets of primitive vari¬ 
ables. We will write the interface state generically as: 

qs = qs-^2(l {v) - AqMyM (Al) 

V 

where s refers to the left or right state at an interface and the sum includes only those waves 
moving toward that interface. We introduce the shorthand for a jump: 

A qM=q,-X<">(q) (A2) 


A.l. Standard primitive variables 


The standard set of primitive variables used in CASTRO for the characteristic tracing 
and prediction of the interface states are q = (p,u,p, pe) J . This evolution is governed by 
Eqs. 0a] to 03 and [HI Written in the form 


<k + A(q)<h = o 


(A3) 


the matrix A takes the form 


A (q ) 


( u p 0 0 ^ 

0 u 1/p 0 

0 pc 2 u 0 

^ 0 ph 0 u j 


(A4) 


with eigenvalues A^) = u — c, A A = u, AA = u, and Al + ) — u + c. We denote the 
second instance of the u eigenvalue with the superscript e, and its associated eigenvector 
only appears when the primitive variable system is augmented with pe. 


Expressing the eigenvectors as a matrix, R 
we have 

/ 1 10 1 \ 


( r ( )| r (°)| r ( e )| r (+)) and L = (/*■ )|/A|/( e )|/(+)jT ; 
/ 0 -p/2c 1/2c 2 0 \ 


—c/p 0 0 c/p 

c 2 0 0 c 2 

\ h 0 1 h j 


L 


1 0 -1/c 2 0 

0 0 -h/c 2 1 

\ 0 p/2c 1/2c 2 0 ) 


(AS) 


The derivation of these follows from Ar^ = \(A r { v ) and l^A = AA)/0) ; with /A • r A) = Sij. 
Note that since the u eigenvalue is degenerate, A = A qi e \ 
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We can now write out the form of the interface states takes by simply multiplying out 
the dot products and doing the sum. We introduce the following notation (based on CW): 


/?(-) = ((<->. A 9 <->) = £(-A U (-> + ^kb (A6a) 

/?<”> = (fM • A„W) = Ap (=)_^hl (A6 b) 

cr 

AW = (((')• A«P) = + A(/>e)<”> (A6 c) 

c 

& (+) = (/ (+) • Ag«) = £ f Aw« + (A6d) 


For these /3’s and those that follow in the other sections, if the respective wave is not moving 
toward the interface, then the associated is set to 0. With these definitions, and Eq. IA11 
we have 


Ps 

= (5-(/?<-'+ /3W + /3W) 

(A7a) 

Ug 

= 5 - (--/?<-> + 

VP P J 

(A7b) 

Ps 

= P - (c 2 /?'-' + e 2 ,3< +) ) 

(A7c) 

(pe) fl 

= (73 - (hAF + A< e) + AA< +) ) 

(A7d) 


A.2. Specific volume in place of density 


CW consider a system with t = 1/p replacing p in the primitive variable system and 
do the characteristic tracing in terms of this. We will denote this system as q— (r,u,p, e) T . 
Now Eq. |4a] is modified by substituting p = r” 1 and expanding: 


And our system is 


dr dr du 
dt U dx dx 


q t + Aq x = 0 


A(q) 


( u —t 0 0^ 

0 U T 0 

0 c 2 /r u 0 
^ 0 pr 0 u j 


(A8) 

(A9) 


(A10) 


with 
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This has the same characteristic polynomial, |A — A/| = 0 as A(q), so the eigenvalues are 
unchanged. The matrices of eigenvectors are: 


/ 1 1 

c/r 0 
-c 2 /r 2 0 

\ -p o 


0 

0 

0 

1 


1 \ 

—c/r 
—c 2 /r 2 

~P 


L 


( 0 r/2c 
1 0 

0 0 

^ 0 —r/2c 


-r 2 /2c 2 0 \ 
r 2 /c 2 0 

-pr 2 /c 2 1 
- t 2 / 2 c 2 0 ) 


(AH) 


As above, 


we introduce the notation 




(Z°M ■ AgM): 


= (f H ■ A4<-») 

/?h = (! ( °> • Ag<°>) 

'0P = (i (e) • AjW) 

/?i +) = (f <+) ■ A<j<+>) 


1 

2 C 

Arf o) + 


Ani } - 

A pi o) 

C ' 2 


Ap, 


(-)' 


C 


+ Ae<"> 

£<2 s 


1 

2 C 


-Am( +) - 


Ap, 


(+)’ 


C 


(A12a) 

(Al2b) 

(Al2c) 

(Al2d) 


These /3’s (excluding the ‘e’ case) are identical to Eq. 3.7 in CW. Inserting into Eq. IA1I 
gives: 


T. = f,-(/§<-> + /§<°> +,3< +) ) (A13a) 

u, = - (c/3<-> - C/§i+>) (A 13b) 

P, = P, - (-C 2 /3‘-> - C 2 /3‘+>) (A13c) 

e. = e. - (-jbh + - P/3i +l ) (A13d) 


If we express r = 1/p and Is°\t) = 1/X ( f\p), then Eqs. IA13al to IA13cl are identical to 
CW Eq. 3.6, giving: 

-1- (s'-'+Z^+zW)] (A14) 


Ps = 


This substitution says that we do the parabolic reconstruction of p, but the characteristic 
projection with r. We note that in evaluating these expressions, CW use C = \JTiftp instead 
of C constructed from the time-level n data. 
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A.3. Using 7 e 


For a general equation of state, we can combine the ideas of CG and CW to predict 7 e 
on the interfaces by first reconstructing it as a parabola. This is in spirit of what the FLASH 
dimensionally-split PPM solver does, but as we’ll show here, there is a correction term that 
is necessary to account for the jump in y e across all waves. 


We consider the same system as CW, augmented with the evolution equation for 7 e . 
We’ll consider the primitive variable system q = (r, u,p, 7 e ) T . Using Eq. [TJ our system can 
be written as: 


q t + Aq x = 0 


(A 15) 


with 


A(q) = 


( u 
0 
0 

V° 


—T 

U 

c 2 /r 

—a 


0 

r 

u 

0 


0 \ 
0 
0 


u 


(A16) 




where we write a = (% — l)( 7 e —Ti) for shorthand. We notice that the structure and elements 
of this matrix are identical to that of A(q ) with the change pr —>■ —a. The eigenvectors are 
then easily computable and found as: 


/ 

1 

1 

0 

1 

\ 


( 

0 

r/2c 

—r 2 /2c 2 

0 

\ 


c/r 

0 

0 

—c/r 




1 

0 

r 2 /c 2 

0 






L = 







—c 2 /r 2 

0 

0 

-c 2 /r 2 




0 

0 

«r/c 2 

1 


V 

a/r 

0 

1 

a/r 

) 


V 

0 

—r/2c 

—r 2 /2c 2 

0 

/ 


(A17) 


As above, we introduce the notation = (7^ • A qA). Because of the similarities to 
the CW system derived above, we see that /3s~ ) = (3s~\ As 0) = (3s°\ and ^i +) = As + \ and 

= ((<«> • A$<«>) = + A 7 «M (A18) 

Considering only how the q e state appears in the characteristic tracing, we see: 

(A19) 

r r 

and if we simplify things by taking the reference state, 7 es to be zero, then we find: 

/3 ( “) Ji 0) (p) ft + ) 

+ tC 2 


les = A 0 \le) + a 


r 


r 


(A20) 
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We note that with a constant-gamma gas (like an ideal gas), then zi° ( j e ) = y e and a = 0, as 
expected. Written as above, we see that the % interface state can be viewed as the average 
of 7 e under the reconstructed parabola over the region traced by the = u wave plus a 
correction term that is proportional to a. In FLASH, to the best of our knowledge, only the 
first term is present in the dimensionally-split PPM solver. 


B. Exact Riemann Solution with a Degenerate Gas 

To test the variations on the standard PPM method we need a test problem with an 
exact solution for the general stellar equation of state. CG (section 1) describe how to exactly 
solve the Riemann problem for a general equation of state. As with a gamma-law gas, across 
the left and right waves, different functions connect the left/right state to the star state, 
and the resulting equation for p* needs to be solved via an iterative procedure. However, 
unlike the gamma-law gas, one cannot write down a closed-form expression for either the 
rarefaction or shock cases. For the rarefaction, a system of ODEs must be integrated, while 
for the shock, a root-finding procedure operates on the Rankine-Hugoniot jump conditions. 

Note that in the following, we assume a constant composition. Adding a composition 
jump at the interface is straightforward, as the composition will only change across the 
contact wave (this follows from the structure of the right eigenvectors). Below we summarize 
the solution procedure from CG, filling in some implementation details: 


Initi al gue ss: We need an initial guess for p*. We use the two-shock approximation 
from Toro ( 1997 1. Eq. 9.42: 


P * = [( W rPl + Wip r ) + WtWriut ~ U r )}/(W l + W r ) 


(Bl) 


with the initial guess for the wave speeds taken to be the Lagrangian sound speed, 
W s = V^lsPsps 

• p* iteration loop: We use perform the Newton iteration described in CG. For s = {l, r}, 
we apply the shock jump conditions if p* > p s , and use the Riemann invariants for a 
rarefaction otherwise. In each case, our goal is to find Z s = |dp*/dn* ]S | and W s . 


— Shock solution: For a shock, we solve the Rankine-Hugoniot jump conditions. 
This takes the form (see CG Eq. 12): 


Ws) 


w s 2 



P*-Ps 

w s 2 



\(pl-p 2 s) (B2) 
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with the solution corresponding to f(W s ) = 0. Here, we used the mass Rankine- 
Hugoniot condition to express the density in the EOS function in terms of p* and 
W s : 


P* 


J_ _ p* -p 8 
,P. w? . 


(B3) 


We solve this via Newton iteration, using the derivative: 


f(W s ) = 2W s {e(p*,P*)~e s } 


2 de 

W~ s Tp 


(p* 


Ps)pI 


(B4) 


We iterate until \f/f'\ < eW s , where e is a small tolerance. We evaluate that 
thermodynamic derivative by expressing our EOS as e(p, T(p, p)): 


de 

de 

de 

dT 

de 

de 

dp 

f dp 


dp 

P d P 

T dT 

P d P 

P d P 

T dT 

P d P 

t \ dT 

J 


(B5) 


where we used the constancy of pressure as: dp = (dp/dp)\Tdp + (dp/dT)\ p dT = 0 
to eliminate dT/dp | p . 

Once we End W s , we can evaluate Z s using CG Eqs. 20 and 23. Two more 
thermodynamic derivatives are needed: 


dp 

de 

dp 

dp 


dp 

f de 

V 

~ dT 

l &T 

j 

p 

p \ 

pj 

dp 

de 

( 

e d P 

t dp 



-i 



and in terms of specific volume, r = 1 /p, 


dp 


T 



(B 6 ) 

(B7) 


(B 8 ) 


These can be derived in a similar fashion as shown above. 

We use the approximation for W s in the case where we know 7 * (CG, Eq. 24) as 
the initial guess for W s . We estimate a value for 7 * from the evolution equation 
CG derive, CG Eq. 31. 

— Rarefaction solution: In the case of a rarefaction, we simply integrate the Riemann 
invariant ODEs. We integate: 


dr 

\-c - 2 

left rarefaction 

dp 

\-G - 2 

right rarefaction 


du 

dp 


—C 1 left rarefaction 
C ' 1 right ratefaction 


(B9) 
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Here C = C(r,p) is the Lagrangian sound speed. The integration is done simply 
with 4-th order Rnnge Kutta from p = p s to p*. Then W s = C and Z s is given 
by CG Eq. 16. 


• Sampling the solution: 

Once we have solved for the wave structure, we can sample the solution at a given 
time, t. Given N points, the coordinate centers are ay = (i + 1/2) Ax for % — 1,..., A, 
with Ax = (x max — x min )/N. Then we define 


6 


ay 


U'max 


(BIO) 


The speed of the contact n* tells us which states we need to deal with, and the evalu¬ 
ation of the solution proceeds from CG Eq. 15 and following. 

The only complication is the case where we are sampling inside the rarefaction. As 
suggested in CG, we integrate the Riemann invariants, however, we found that it is 
easiest to rewrite them with u as the independent variable. We integrate: 


dr j C 1 left rarefaction dp 

du 1 —C~ l right rarefaction du 


—C left rarefaction 
C right rarefaction 


We integrate from u = u s to 


(BH) 


^stop 


£ + c(r,p) 
f-c(r,p) 


left rarefaction 
right rarefaction 


(B12) 


We don’t know u stop at the start of the integration (since c depends on the state), so 
we evaluate it each step. Again we use Runge-Kutta integration, and we adjust our 
stepsize to stop at the converged u s t op - 


We use this solution to generate several exact shock tube profiles to compare with the 
PPM solutions. It is easier to specify the initial conditions in terms of a jump in ( p,u,T) 
instead of (p,u,p) and then compute the pressure via the equation of state. Additionally, 
for this equation of state, we need to specify the composition of the gas—we take it to be 
pure 12 C. 

Since even different versions of the same equation of state can have slightly different 
thermodynamics, we do not provide a table of the solution output for the problems listed 
here. Instead, the complete source code for this exact solver is available in the CASTRO 
source distribution (in the Util/exact_riemann/ subdirectory). This can be built with 
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either the general stellar equation of state flTinnn es fe Swestvl 2000) or a gamma-law EOS. 
All the necessary inputs hies for the tests described here are provided. 


We devise 4 tests. The first 3 are analogous to tests 1 through 3 in Toro (119971 )): test 1 is 
a Sod-like problem, test 2 is a double rarefaction, and test 3 is a strong shock. The last test, 
test 4, mimics the conditions present at the edge of a star when it is not very well resolved 
on a grid—a common issue with multi-dimensional stellar simulations. Table Q] lists all the 
parameters, including the length of the domain, L , and the end time of the simulation, t. 
All tests used outflow (zero-gradient) boundary conditions. For test 4, because of the large 
drop in density, the Coulomb corrections in the equation of state were disabled to avoid 
unphysical regions in the thermodynamic space. 
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Table 1. Test problems. 


quantity test 1 test 2 test 3 test 4 



left 

right 

left 

right 

left 

right 

left 

right 

p (g cm 3 ) 

10 7 

10 6 

10 7 

10 7 

CD 

O 

t— H 

10 6 

10 2 

10" 4 

T (K) 

10 8 

10 6 

10 8 

o 

00 

o 

CD 

10 6 

10 7 

10 7 

u (cm/s) 

0 

0 

1 

o 

00 

O 

00 

0 

0 

0 

0 

L (cm) 

10 6 

10 5 

2 x 10 5 

10 5 

t (s) 

8 x 10" 4 

8 x 10” 5 

2 x 10" 4 

3 x 10" 4 






pressure (erg/cc) 
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Fig. 3.— Sod-like test 1 problem with the general EOS comparing the differences between 
the CG and CGF Riemann solvers 
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Fig. 4.— Errors in the Sod-like test 1 problem with the general EOS comparing the differ¬ 
ences between the CG and CGF Riemann solvers 
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x xlO 5 

Fig. 5.— Solutions for the “stellar edge” test 4 problem with the general EOS comparing 
the differences between the CG and CGF Riemann solvers 
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xlO 6 


Fig. 6.— Sod-like test 1 problem with the general EOS comparing the use of r = 1/p (CW) 
to p (MC) in the characteristic tracing. Only the temperature is shown 
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Fig. 7.— A comparison of CASTRO to FLASH for test 1. 
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x xlO 5 X xlO 5 



Fig. 8.— A comparison of CASTRO to FLASH for test 2. 
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Fig. 9.— A comparison of CASTRO to FLASH to test 3. 
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x xlO 5 x xlO 5 



Fig. 10.— A comparison of CASTRO to FLASH to test 4. 






































